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We study the effects of nonzero time delays in stochastic synchronization problems with linear 
couplings in complex networks. We consider two types of time delays: transmission delays between 
interacting nodes and local delays at each node (due to processing, cognitive, or execution delays). 
By investigating the underlying fluctuations for several delay schemes, we obtain the synchroniz- 
ability threshold (phase boundary) and the scaling behavior of the width of the synchronization 
landscape, in some cases for arbitrary networks and in others for specific weighted networks. Nu- 
merical computations allow the behavior of these networks to be explored when direct analytical 
results are not available. We comment on the implications of these findings for simple locally or 
£f} . globally weighted network couplings and possible trade-offs present in such systems. 
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I. INTRODUCTION 



Since the classic works by Kalecki [l] and Frisch & Holme Q on the emergence of macro-economical patterns 
(business and economics cycles), it has been well-established that time delays occurring on microscopic scales can 
have profound effects on the global response of complex systems. Among other early key results were the works by 
Hutchinson [3} and May Q, showing that time delays can have fundamental impact on logistic growth in population 
dynamics jEfl. The importance of time delay becomes even more explicit in interacting individual- or agent-based 
models |6l4l0l|. where the delays can correspond to time scales in the interactions (e.g. transmission delays) or to 
time scales of the local decision and execution by the individuals. In this paper, we consider the simplest - yet 
fundamental - model for such networked systems, taking into consideration the effects of the network topology and 
couplings, noise, and time delays [ill ]. This paper provides an extended account of our recent Letter [il| . providing 
more details, generalizations and comparisons to certain weighted networks, and considering different types of delays. 

In network synchronization [l2j . coordination, or consensus problems individuals or entities represented by 
nodes in the network attempt to adjust their local state variables (e.g., pace, load, phase, or orientation) in a 
decentralized fashion. (In this paper, we use the terms synchronization, coordination, and consensus synonymously 

■ in this broader sense.) Nodes interact or communicate only with their local neighbors in the network, often with the 
, intention to improve global performance. These couplings can be represented by directed or undirected, weighted or 

unweighted links. Applications of the corresponding models range fro m p hysics, biology, computer science to control 
theory, including synchronization pro blems in distributed computing [131 ] , symbolic dynamics |14j , congestion control 
in communication networks [E [l5l4l8j | and in vehicular traffic pE H3 > flocking animals [2l| - |23j , bursting neurons [13] , 
' and cooperative control of vehicle formation pBj . 

■ Synchronization, coordination, or consensus in complex networks cuts across numerous fields that address global 
behavior through decentralized local actions facilitated by sparse interactions. There has already been much in- 
vestigation into the efficiency and optimization of synchronization [l2l l2fj| - [29j in weig hted HE EE 13 and directed 
[E III III topolo gies. Because of limitations in communication, transportation, processing, or cognitive resources, the 
local information on the state of the network neighborhood may not always be current, nor is it even given for the 
same instant at a past time for all components. These time delays can have drastic effects on system behavior Q and 
further complicate predictability of the network's global performance. 

The impact of time delays on stochastic differential equations involving a single stochastic variable, with recent 
applications to postural sway [33l |34| , stick balancing at a fingertip [35l . [3rl . and the scaling of congestion window 
in internet protocols [37]], have been investigated in the past two decades 38 41]. Here, we focus on the interplay 
of network topology, couplings, noise, and time delays. Our motivation is to understand how network-connected 
individuals contribute to global goals by performing delayed actions and/or using delayed information facilitated by 
local interactions in a noisy environment. 
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The phenomena of spontaneous synchronization, coordination, or consensus arise in a variety of disciplines 0, @, 
[l6j |. For example, it describes the consensus that arises in bird flocks as each bird makes velocity adjustments to match 
the group, which is crucial in accomplishing such tasks as avoiding predators [2 2| . [2 3l|. Similarly, it can be applied 
to a collection of autonomous vehicles working cooperatively to carry out a task |l8[ . Risk can be managed without 
central governance in uncertain environments through the synchronicity or spontaneous cooperation of individuals. 
This appears in economics when considering stock trades [421 ] ; in ecology there is the reward of reproduction and the 
danger of predation for chirping cicadas and flashing fireflies [43| . While there are adversarial relationships between 
individual participants, there are still mutual benefits (predictive insight or bodily protection) from the collective 
behavior. Massively parallel and distributed computing schemes require synchronization across processors [l3|, l44T - |46j 
in order to avoid diverging progressions of simulation time but must be balanced with the cost of communication. 
Synchronization of coupled phase oscillators [13] (the Kuramoto model (48|) has many applications, recently to spatial 
patterns in flashing microfiuidic arrays [491 ] and to circuits comprised of optomechanical arrays [50|. In neural networks, 
time delays critically affect the synchronization of excitatory fronts [5ll453| . All these examples are instances of a 
group coming to consensus [(| without an omniscient global operator. They fundamentally rely on the communication 
between individuals, which may be (and often is) sent through noisy channels [39 41]. 



A. The Model 



In the model we consider here, the state of each node i is described by a local scalar state variable hi- In stochastic 
network coordination/consensus problems, nodes locally adjust their state in an attempt to match that of their 
neighbors through linear couplings in the presence of noise. However, they react to the information or signal received 
from their neighbors with some time lag, and the evolution of the states of the nodes is governed by the differential 
time-delay equations 

dMt) = - ]T CaMt - t?) - hj(t - r? - r£)] + m(t) . (1) 

i 

Here, Cij is the coupling strength between nodes i and j, and rji is the noise present at node i, satisfying (T]i(t)r]j(t')) — 
2DSijS(t — t'), where D is the noise intensity. In general, the time delays can be heterogeneous, depending on the 
properties and network locations of both nodes: r° is the local delay at node i, corresponding to processing, cognitive, 
or execution delays, while r« is the transmission delay between nodes i and j. Without the noise term, the above 
equation is often referred to as the (deterministic) consensus problem 0, [l8| on the respective network. In this sense, 
the networked agents try to coordinate or reach an agreement or balance regarding a certain quantity of interest. 
A standard measure of synchronization, coordination, or consensus in a noisy environment is the width [T3 . HH 

(« 2 (*)> = (iE[M*)-^(*)] a ). ( 2 ) 

\ i=l ' 

where h(t) = (1/N) J^Zi MO is the global average of the local state variables and (. . .) denotes an ensemble average 
over the noise. A network is "synchronizable" if it asymptotically reaches a steady state with a finite width, i.e. 
(w(oo)) < oo. When the network is well synchronized (or coordinated), the values hi for all nodes are near the global 
mean h and the width is small. 



B. Coordination without Time Delays 

Without time delays, Eq. (fT]) takes the form 

dMt) = -Y, C v\ h i{t) - hj(t)]+in(t) = -J2 T ijhj(t) + (3) 

3 3 

where T^j = 5y Cu — Cij is the network Laplacian. Eq. (J3|) is a multivariate Ornstein-Uhlenbeck process [54{ and 
is also referred to as the Edwards- Wilkinson process [55[ on a network [l3|> [l5[ . Starting from a flat initial profile 
{hi(0) — 0}^Li for symmetric couplings, one can show that the width evolves as [H| 

k=l " 



2 



where Afc, k = 0, 1, 2, . . . , N — 1, are the eigenvalues of the network Laplacian. Note that as a result of measuring the 
local state variables hi from the mean h in Eq. ([2]), the singular contribution of Ao = (associated with the uniform 
mode) automatically cancels out from the sum in Eq. ((4]). Thus, a finite connected network is always synchronizable 
with steady-state width 

<»"«»>»- < 5 > 

fc=i 

In the limit of infinite network size, however, network ensembles with a vanishing (Laplacian) spectr al g ap may 
become unsynchronizable, depending on the details of the small- A behavior of the density of eigenvalues 12t fl3l Il5j. 



This type of singularity is common in purely spatial networks (in particular, in low dimensions) where the relevant 
resp onse functions and fluctuations diverge in the long- wavelength (small- A) limit [H, [56|. In complex networks 
[57H60( these singularities are typically suppressed as a result of sufficient amount of randomness in the connectivity 
pattern generating a gap or "pseudo" gap. (13, 12111 [olTI&ij ■ 

As is also clear from Eq. ([5]). synchronization or coordination can be arbitrarily improved in this case of no time 
delays, e.g., by uniformly increasing the coupling strength by a factor of a > 1, resulting in Gy —¥ aCij (A& — ¥ aXk) 
and yielding 

(w 2 (oo)) a = I(™ 2 (oo)) CT=1 . (6) 
er 

The stronger the effective coupling a (e.g., achieved by more frequent communications in real networks), the better 
the synchronization; the width is a monotonically decreasing function of a. 



II. UNIFORM LOCAL TIME DELAYS 



We first consider the case with symmetric coupling Cy = Cji when transmission delays are negligible (tjJ = 0) and 
local delays are uniform (r° = r). Then Eq. (jTJ is governed by a single uniform time delay 



N N 



dMt) = -J2 Cij[hi(t - r) - hj(t - t)] + m(t) = -J2 r 'A (* - T ) + m(t) ■ (7) 

This equation has a similar form to that of Eq. ([3]) but with the inclusion a delay r. 



A. Eigenmode Decomposition and Scaling 

By diagonalizing the symmetric network Laplacian T, the above set of equations of motion decouples into separate 
modes 

d t h k (t) = -X k hk(t-T)+f} k {t) , (8) 

where A& (k = 0, 1, 2, . . . , N — 1) are the eigenvalues of the network Laplacian, and and fjk are the time-dependent 
components of the state and noise vectors, respectively, along the k-ih eigenvector. Thus, the amplitude /i& of each 
mode (with the exception of the uniform mode with Ao — 0) is governed by the same type of stochastic delay-differential 
equation 

d t h(t) = -Xh(t - r) + fj(t) , (9) 

with A > 0, where we temporarily drop the index k of the specific eigenmode for transparency and to streamline 
notation. 

While the above stochastic delay-differential equation has an exact stationary solution for the stationary-state 
variance [38|, HTJ] , we first review the formal solution [TJ |(55[ which provides some insights and connections between 
the solutions of the underlying characteristic equation and the existence (and the scaling) of the stationary-state 
fluctuations of the stochastic problem. The formal solution can also be applied to more general linear (or linearized) 
coordination problems with multiple time delays [66| . and can serve as the starting point to extract the asymptotic 
behavior [67| near the singular points (synchronization boundary). 
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FIG. 1: (Color online) Time evolution of an individual mode obtained by numerically integrating Eq. Q with A=l, D=l, and 
At=0.001 for several delays chosen to show the various behaviors across the separating/critical points Ar=l/e and n/2; (a) 
Ar=0.2 < 1/e, (b) 1/e < Ar=1.5 < tt/2, and (c) Ar=1.7 > tt/2. 



Performing standard Laplace transform on Eq. © [with /i(s) = / Q °° e st ft,(i)c?f], the characteristic equation associ- 
ated with its homogeneous (deterministic) part becomes 

g(s) = s + Xe ST = . (10) 

As shown in Appendix |A1 (with h(t) = for t < 0) the time-dependent fluctuations can be written formally as 

-2D(1 - e^+^fj)*) 



a,/3 



g'(sa)g'(sp)(s a + Sp) 



Hence, they remain finite (i.e., a stationary distribution exists) if 

Re(s Q ) < , 



(11) 



(12) 



for all a, where s a , a = 1, 2, . . ., are the solutions of the characteristic equation, Eq. (|10[) . on the complex plane. We 
can explicitly make the simplification 



-2D 



-2D 



n g'{s a )g'{sp)(s a + Sf}) (l + TS a )(l + TS )(s a + Sp) 



(13) 



Eq. (flU|) is perhaps the oldest and most well-known (transcendental) characteristic equation from the theory of delay- 
differential equations @, ID, [TH, |68| , with the linear stability analysis of numerous nonlinear systems reducing to this 
one. It has an infinite number of (in general, complex) solutions for r > and the condition in Eq. (1121) holds if 



At < tt/2 



(14) 



Long-time dynamics of the solution of Eq. (j9]) is governed by the zero(s) of Eq. (fl0|) with the largest real part. In 
particular, for At < 1/e, the zero with the largest real part is purely real, hence no sustained oscillations occur 
[Fig. (Ha)]. For 1/e < At < tt/2, all zeros have imaginary parts (including the ones with the largest real part) and 
are arranged symmetrically about the real axis. This results in persistent oscillations that do not diverge so long as 
condition (TT4l is satisfied, as shown in Fig.QJb). The first pair of zeros to acquire positive real parts are the two with 
smallest imaginary parts. Once the product At fails to satisfy the condition in Eq. (I14[) , the oscillation amplitude 
grows in time [Fig. HJc)]. Specific time series for (h 2 (t)} are shown in Fig. [21 where the real parts of solutions have 
become positive for delays t = 1.60 and 2.00 but remain negative for the rest. 
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FIG. 2: (Color online) Time series of the fluctuations of a single mode (A=l) averaged over 10 4 realizations of noise (with 
D=l) by numerically integrating Eq. ((9} with A£=0.01 for different delays (from bottom to top in increasing order of r). 



To obtain the general scaling form of the fluctuations in the stationary state, we define z a = rs a (a — 1,2,...). 
One can easily see that the new variables z a are the corresponding solutions of the scaled characteristic equation, 



z + Are 2 = , 

and hence can only depend on Ar, 1.0. — 

(At). Thus, 



Substituting this into Eq. (fT3]) yields 



where 



s Q (A,r) = -z a (Xr) . 
r 



(h 2 (oc)) = D T f(\r) , 



-2 



a,/3 



(1 + Z a )(l + Zp){z a + Zfj) 



(15) 



(16) 



(17) 



(18) 



is the scaling function. This scaling [Eq. (|T7| ] is illustrated by plotting (h 2 (oo)) /r vs Ar, fully collapsing the data for 
different r values (with fixed noise intensity D) [Fig. [3]. 

As mentioned earlier, Eq. Q has an exact solution for the stationary-state variance obtained by Kiichler and 
Mensch [381 ] (briefly reviewed in Appendix [Bj , providing an exact form for the scaling function 



/(Ar) 



1 + sin(Ar) 
Ar cos(Ar) 



(19) 



The asymptotic behavior of the scaling function near the singular points, Ar — and Ar = n/2, can be immediately 
extracted from the exact solution given by Eq. (|19p (see also Ref. [67| for a more generalizable method) , 



/(Ar) 



1 

At" 



7r(7r/2 - Ar) 



< Ar < 1 



< - - Ar < 1 . 



(20) 



The scaling function f(x) (x = Ar) is clearly non-monotonic; it exhibits a single minimum, at approximately x* 
0.739 with /* = f(x*) « 3.06, found through numerical minimization of Eq. (fT9"|) . The immediate message of the 
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FIG. 3: (Color online) (a) Steady-state fluctuations of an individual mode as a function of A obtained by numerical integration 
of Eq. for several delays with D=l and A£=0.01. (b) Scaled fluctuations of an individual mode and the analytic scaling 
function Eq. (|19l) (solid curve). 



above result is rather interesting: For a single stochastic variable governed by Eq. @ with a nonzero delay, there is 
an optimal value of the "relaxation" coefficient, A* = x*/t, at which point the stationary-state fluctuations attain 
their minimum value (/i 2 (oo)) = Drf* ss 3.06-Dr. This is in stark contrast with the zero-delay case (the standard 
Ornstein-Uhlenbeck process [HJ]) where (h 2 (oo)) = D/X, i.e., the stationary-state fluctuation is a monotonically 
decreasing function of the relaxation coefficient. 

B. Implications for Coordination in Unweighted Networks 

Since the eigenvectors of the Laplacian are orthogonal for symmetric couplings, the width can be expressed as the 
sum of the fluctuations for all non-uniform modes 

< w2 m> = ^ E = -jf E /( A * T ) > ( 21 ) 

k=l k=l 

where Afe is the eigenvalue of the fcth mode. Thus, condition (fl"4|) must be satisfied for every k > mode for 
synchronizability, or equivalently (69[, 

7T 

X m ax T < ~ • (22) 

The above exact delay threshold for synchronizability has some profound consequences for unweighted networks. Here, 
the coupling matrix is identical to the adjacency matrix, Cy = Ay, and the bounds and the scaling properties of the 
extreme eigenvalues of the network Laplacian are well known. In particular [70l [7l| , 

N 



N - 1 



where /c max is the maximum node degree in the network [i.e., (A max ) = 0((fc max ))]. Thus, rfc max < 7r/4 is sufficient for 
synchronizibility [69[, while rfc max > n/2 leads to the breakdown of synchronization with certainty. These inequalities 
imply that even a single (outlier) node with a sufficiently large degree can destroy synchronization or coordination in 
unweighted networks (regardless of the general trend, if any, of the tail of the degree distribution). Naturally, network 
realizations selected from an ensemble of random graphs with a power-law tailed degree distribution typically have 
large h ubs, making them rather vulnerable to intrinsic network delays 0, HH. For example, Barabasi- Albert (BA) 
[58|,[59[ and uncorrelated [72|,[23| scale-free (SF) networks with structural degree cut-off (yielding A max ~ fc max ~ N 1 / 2 ) 
and similarly, SF network ensembles with natural cut-off (exhibting A max ~ & max ~ 

jVV(7-i)) for N > j [63,(71), 

are particularly vulnerable. Thus, for any fixed delay, increasing the size of scale free networks will eventually lead 
to the violation of condition (|22[) . and in turn, to the breakdown of synchronization. In contrast, the typical largest 
degree (hence the largest eigenvalue of the Laplacian) grows much slower in Erdos-Renyi (ER) random graphs [74| . 
Ln(iV). 
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FIG. 4: (Color online) The fraction of synchronizable networks p s (r, N) taken from ensembles of 10 4 random constructions of 
ER and BA networks with (k) « 6. (a) p s vs. N. (b) and (c) are scaled plots of the same data according to Eq. (|24l) . for ER 
and BA networks, respectively. 



To illustrate the above finite-size dependence, we define the fraction of synchronizable networks p s (r, TV), which is 
equivalent to the probability that a randomly chosen realization of a network ensemble satisfies A max <7r/2r. Thus, 
p s (r,N) = P^(it/2t), where P^(x) is the cumulative probability distribution of the largest eigenvalue of the network 
Laplacian. In Fig. 21 we show the fraction of synchronizable networks for BA and ER network ensembles by employing 
direct numerical diagonalization of the corresponding network Laplacians and evaluating condition (|22[) for each 
realization. For N 1 the cumulative distribution for the largest eigenvalue exhibits the asymptotic scaling P^ (x) ~ 
4>{x / (X max (N)}) (64[. Thus, the fraction of synchronizable networks should scale as 

Ps (t,N) = P$(tc/2t) ~ <^(7r/2r<A max (iV))) = V(r(A max (TV))) . (24) 

In Figs. S](b) and (c), we demonstrate the above scaling for ER and BA networks, respectively. 

Since the scaling function is known exactly [Eq. (|T9l)]. the eigenmode decomposition [Eq. (|21"1) ] allows one to evaluate 
the stationary width for an arbitrary network with a single uniform time delay by utilizing numerical diagonalization 
of the network Laplacian 

/ 2/ sv Dt V^jv\ n £»t ^ 1 + sin(A fc r) 

(w (oo) = — > /(A fc r = — > — - . (25 

The optimal (minimal) width occurs when all eigenvalues of the Laplacian are degenerate so that the couplings and/or 
delay can be tuned to the minimum of Eq. (1191) . For each mode in Eq. (|25[) . such degeneracy is present in the case 
of a fully-connected network with uniform couplings, optimized to Cy = x* /Nt (i ^ j) and Cu — 0. For general 
networks, better synchronization can be achieved when the eigenvalue spectrum is narrow relative to the range of 
synchronizability so that most eigenvalues can fall near the minimum of Eq. (|19p . We have not investigated in detail 
how to achieve a narrow spectrum, but strategies for doing so by tuning coupling strengths or adding/removing links 
have been explored by others [HI, HH . 



C. Scaling, Optimization, and Trade-offs in Networks with Uniform Delays 

With the knowledge of the scaling function in Eq. (|2"5)l , one also immediately obtains the width for the case of an 
arbitrary but uniform effective coupling strength a, where Cy = cAij. The effective coupling strength can now be 
tuned for optimal synchronization. However, there is a trade-off between how well the network synchronizes and the 
range over which it is synchronizable. When the eigenvalue spectum is not narrow, diminishing the couplings uniformly 
in order to satisfy Eq. (|2"2"]l may cause small eigenvalues to be pushed farther up the left divergence of the scaling 
function. Figure EJa) shows this trade-off in uniform reweighting (Cy — > aCij). The monotonicity of these widths 
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FIG. 5: (Color online) Stationary-state widths obtained through numerical diagonalization and utilizing Eq. (|25|) for a typical 
BA network with N = 100 (a) for several coupling strengths, (b) for several delays, and (c) scaled so that the nonzero delay 
curves collapse. 



means that the uniform delay should always be minimized to obtain the best synchronization. The same conclusion 
can be drawn from Fig. [5jb) , which shows that networks synchronize better and do not become unsynchronizable until 
greater link strengths when the delay r is minimized. Because globally reweighting the coupling strengths corresponds 
to a uniform scaling of the eigenvalues, we can define the width of a network by a scaling function F(ctt) (see Fig.(5Jc)) 

N-l 

( W 2 M) CT . r = E /( ffA * T ) = DrF(ar). (26) 

Fluctuations from small eigenvalues dominate other contributions to the width for small err, hence the optimal value 
occurs near the end of the synchronizable region, where the network fails to meet condition (|22[) . 

As an alternative to varying the (effective) uniform coupling strength a, consider a scenario where the frequency 
(or rate) of communication is controlled for each node according to 

N 

d t hi(t) = - Pi (t) J2 Aij[h{t - r) - h 3 (t - t)} + m(t) . (27) 

3=1 

In the above scheme, Pi(t) is a binary stochastic variable for each node, such that at each discretized time step, 
Pi (t) = 1 with probability p and pi (t) = with probability 1 — p (for simplicity, we employ uniform communication 
rates). The local network neighborhood remains fixed, while nodes communicate with their neighbors only at rate p 
at each time step. As an application for trade-off, consider a system governed by the above equations and stressed by 
large delays, where local pairwise communications at ratep=l would yield unsynchronizability, i.e., T\ max >n /2 (see 
Fig. [5J)- The width diverges for one of two reasons: either communication is too frequent and the system fails to satisfy 
condition (1221) . or there is no synchronization (p=0) and the system is overcome by noise. However, the divergence 
of the width is faster in the former, accelerated by overcorrections made by each node due to the delay. With an 
appropriate reduction in the communication rate, the width reaches a finite steady state, recovering synchronizability, 
as can be seen in Fig. [5] Decreasing the frequency of communication can counter-intuitively allow a network to become 
synchronizable for delays and couplings that would otherwise cause the width to diverge. 
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FIG. 6: (Color online) Time evolution of the width obtained by numerically integrating Eq. (|27[) with D=l, At=0.005, and 
averaged over 10 3 realizations of noise for several communication rates p on a BA network of size iV=100 and average degree 
(fc)=6 with T\ ma x=l-2 x 7r/2. 

D. Coordination and Scaling in Weighted Networks 

For the case of uniform delays, we compare two cases: networks with weights that have been normalized locally by 
node degree and networks with weights that are globally uniform. The couplings for local weighting are defined as 
Cij = crAij /hi (a common weighting scheme in generalized synchronization problems [HI), while for uniform couplings 
Cij = aAij/(k). In turn, the weighted (or normalized) Laplacian becomes T = <jK^ 1 L where K is the diagonal matrix 
with node degrees on its diagonal, Kij = Sijki, and L is the graph Laplacian, Lij = Sij An — Aij — Sijki — Aij. 
Similarly, for uniform couplings, the corresponding Laplacian becomes T = a(k)^ 1 L. Note that the overall coupling 
strength (communication cost) is the same in both cases, o'Ylij Aij /hi — &J2ij Aij/(k) = <J ^- 

In the locally-weighted case, the eigenvalue spectrum of K~ X L is known to be confined within the interval [0, 2] 
[75j . so any network of this class will be synchronizable, provided ot < 7r/4. With globally uniform weighting, the 
increase of X m ax with N will lead to fewer synchronizable networks as N grows (holding (fc) constant). Figure [TJa) 
shows that it is more likely for an ER network to be synchronizable than a BA network of the same size N when 
the couplings are weighted uniformly by (k) (with all ER networks remaining synchronizable over the range of N for 
the two smallest delays). However, this is not always the case when couplings are weighted locally by node degree 
(Fig. 0(b) ) , although nearly all of these networks remain synchronizable over the delays in Fig. [T^a). The behavior of 
the width for typical networks is shown in Fig. [8] to compare the effects of these two normalizations. In both the BA 
and ER case, synchronization is better and is maintained for longer delays when the coupling strengths are weighted 
locally by node degree. 

III. MULTIPLE TIME DELAYS 

To generalize the basic model, we now allow for a distinction in transmission and processing time delays. In this 
case, Eq. JT]) becomes 

d t hi(t) =~J2 C vMt - T ) - h 3 {t -T - T tr )] + lfc(t) (28) 
3 

where the local delay r Q and the transmission delay r tr are the same for all nodes and links, respectively. Although 
the synchronizability condition and steady state width cannot be determined in a closed form for arbitrary networks 



9 




FIG. 7: (Color online) Fraction of synchronizable networks for (a) uniform global weights and (b) local weights for the same 
ensemble of networks used in Fig. [4] 
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FIG. 8: (Color online) Scaled widths simulated with A£=0.01 of a typical BA and a typical ER network, each of size iV=100 
and with (fc)=6. 



as is the case of Eq. ([7]), focusing on special cases does offer insight. 



10 



A. Fully-Connected Networks 



Consider the case of a fully-connected network of size N with uniform link strengths a, where the local state 
variables evolve according to 

dMt) : 



a 


N- 


1 


a 




N- 


1 


a 




N- 


1 



^2[hi(t - r) - hj(t - t)] + crh t (t - t) - ah^t - 7 t) + ^(i) 

Fijhjit - r) + ahi(t - r) - ah t {t - 7 r) + (29) 

where r = r + rt r , 7 = r /r and = 5ijN — 1. Normalizing the global coupling with l/(iV — 1) assures that the 
coupling cost per node remains constant and the region of synchronization remains finite in the limit of N —> 00. Using 
the fact that the graph Laplacian of the complete graphs has a single, nonzero eigenvalue N [which is (N— l)-fold 
degenerate] , each non-uniform mode (associated with fluctuations about the mean) obeys 

d t h(t) = -ah(t - tt) - j£jK* ~r) + m . (30) 

As in the case of uniform delays, we perform a Laplace transform on the deterministic part to obtain the characteristic 
polynomial and equation, 

9(s) = s + j^je-" + ae-"<™ = . (31) 



Note that for N = 2, the region of stability/synchronizability can be obtained analytically [66|, and for completeness 
we show it in Fig. [9] in the (t ,t) plane. In this simple case of two coupled nodes, the synchronization boundary is 
monotonic, and the local delay is dominant: There is no singularity (for any finite r tr ) as long as trr < 1/2 66], while 
for any rt r , there is a sufficiently large r Q resulting in the breakdown of synchronization. 

For N > 3, the phase diagram (region of synchronizability) can be obtained numerically by tracking the zeros of the 
characteristic equation Eq. (|3Tj) (i.e., identifying when their real parts switch sign) shown in Fig.[9l Note that keeping 
track of infinitely many complex zeros of the characteristic equations would be an insurmountable task. Instead, in 
order to identify the stability boundary of the system, one only needs to know whether all solutions have negative 
real parts. This test can be done by employing Cauchy's argument principle [76l |77| (see Appendix [C] for details). 
Similar to the N=2 case, the local delay is always dominant, i.e., there are critical values of <jt q above/below which 
the system is unsynchronizable/synchronizablc for any r tr . [These critical values approach 7r/2 as N — > 00, since in 
this case Eqs. (|50|) and ([3~Tj) reduce to the familiar forms of Eqs. © and (flU)) , respectively, with the known analytic 
threshold.] The behavior with the overall delay r = t q + r tr , however, is more subtle: There is a range of t d where 
varying r yields reentrant behavior with alternating synchronizable and unsynchronizable regions (as can be seen by 
considering suitably chosen horizontal cuts for fixed t in Fig. [S]). Thus, in this region (for fixed local delays r D ), 
stabilization of the system can also be achieved by increasing the transmission delays. 

In the special case 7 = 0, the network is always synchronizable for all N and the width can be obtained exactly 
(see Appendix IB 21 b). 



1 D(N-l) a + j^- sinh(ar) 



( w 



M> - ]v L ««-» - \ '^^'V^;,, (32) 

with a = cryl — l/(N — l) 2 , as shown in Fig. [TD] For r = r tr — ► 00, the above expression becomes 

2 D(N-1) 1 

(w (00)) = — ; . 33 

V V ^ N CTyT- l/(iV- l) 2 V ' 
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FIG. 9: (Color online) Phase diagram (synchronization boundary) for fully-connected networks with uniform coupling strength 
a/(N — 1) in the (t ,t) plane. (Here, without loss of generality due to scaling, we used <j=1.) (a) All system sizes (for 7V>300 
they are essentially indistinguishable from the iV=30,000 case at these scales); (b) system sizes iV=30, 300, 3,000, 30,000 in an 
enlarged region for visibility. With the exception of the analytically solvable case of N=2 [6(| , the synchronization boundaries, 
corresponding to stability limits, were obtained from the analysis of the zeros of Eq. (1311) . 

B. Locally Weighted Networks 



Now we consider Eq. (128|) with specific locally weighted couplings (already utilized for uniform local time delays in 
Sec. II. D), dj — aAij/ki. The set of differential equations then have the form 



dMt) 



3 

= ~y Lijh^t - t) + ahi(t - t) - ah^t - jt) + rji(t) 



3 



-n- > Tijhj(t - t) + o7i,(f — r) - ahi(t - jt) + r)i(t) 



3 



(34) 



where a controls the coupling strength and V = K 1 L is now the locally weighted network Laplacian (Kij — Sijki, 
and Lij — #y An — Aij = Sijki — Aij). Diagonalization yields 



d t h k (t) = (7(1 - Xk)h k (t - t) - ah k {t - yr) + fj k (t) 



(35) 



where \ k is the eigenvalue of the A:th mode of the normalized graph Laplacian K~ 1 L. Figure [TT] shows the evolutions 
of a particular mode with delays on either side of the critical delay. The characterisitic equation for the fcth mode is 
then 



-*yrs 



g k (s) = s + cr(A fe - l)e TS 4 
Defining the new scaled variable z = ts, this equation becomes 

z + (<tt) (Afe - l)e~ z + (trr)e-T* = 







(36) 



(37) 



Hence, the solutions of the original characteristic equation depends on a and r in the form of Sfc Q = t z ka (crT). 
Although the scaling function of the width in the case of locally normalized couplings with two time delays cannot be 
expressed in a closed form, the general scaling behavior is identical to Eq. (I26p [as follows from the formal solution 
shown in Appendix IDl Eq. (|D11[) ]. i.e., (w 2 (oo)) a , T — DtF(<tt). The corresponding scaling behavior and scaling 
collapse, obtained from numerical integration of Eq. (|34[) . are shown in Fig. 1121 
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FIG. 10: (Color online) Analytic results for stationary-state widths for fully-connected networks of several sizes for the special 
case 7=0 [Eq. p2])]. Here, D=l and <r=l. 

The stability /synchronization boundary was again determined by employing Cauchy's argument principle [7(1 
r77j . applied separately for each mode (Appendix [C]) . Figure [T31 shows the most important eigenvalues to determine 
synchronizability: the greatest restriction to the critical delay r c = (r Q + T tr ) c for a given 7 belongs to either the 
smallest or largest eigenvalues. An alternative presentation is given in Fig. 1 141 which shows that it is not always 
the same eigenvalue that consistently limits synchronizability for all values of 7; rather it is the eigenvalue that falls 
on the lowest point on the boundary curve. The contributions of a few example modes to the width are shown in 
Fig. ITST a). Note that the order of divergences is not the same as the ordered eigenvalues, in accordance with Fig. [T3] 
The contributions of a single mode for various values of 7 is shown in Fig. [T5lb). Since it is r Q that has a greater 
impact on whether or not a network can synchronize, larger total delays r are tolerated for smaller 7 since more of the 
delay comes from transmission. Because of the great sensitivity of (h 2 ) on At near the divergence for longer delays, 
an adaptive algorithm was implemented, which would halve At until consecutive runs agreed within 1%. 

With this understanding of the underlying modes, let us return to synchronization of the entire system. Incorpo- 
rating all relevant eigenvalues results in the synchronization boundary shown in Fig. I16f a) for several representative 
networks. The cut for a carefully chosen local delay in Fig. ITST b) shows the previously mentioned reentrant behavior 
as the transmission delay is increased. Note that the optimal width within each synchronizable region worsens with 
larger delay, so that while synchronizability can be recovered with increasing Tt r , better synchronization is possible 
by decreasing Tt r - To compare the contribution of modes within the synchronizable regime, consider again the two 
topologies of BA and ER graphs. For fixed 7, Fig. [T71 shows that a BA graph remains synchronizable for larger delays 
than a ER graph when the link strengths are weighted by node degree. However, the ER graph synchronizes slightly 
better for the majority of the time that it is synchronizable. Here it is not the topology but the ratio 7 that has the 
most drastic effect. 

When 7 < 1, the mode corresponding to Aq — includes self-interaction terms and has the critical delay 



While the uniform mode does not contribute to the width because h is removed from the state of the network (see 
Appendix [D]) , a diverging mean can introduce egregious truncation errors into the numerical integration if h diverges 
exponentially while the width remains finite. Fortunately, this can be avoided by simulating the network in the 
subspace lacking the zero mode by removing the mean from each time slice. Since the uniform mode is not allowed 
propagate, it does not cause any problem with finite precision. The locations of the zeros' real parts for Eq. (1371) are 
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FIG. 11: (Color online) Time series of the fluctuations of a single mode for several delays obtained from numerical integration 
of Eq. (|35[) with 7=0.5, A=1.8, D=l, and At=0.01, averaged over 10 3 realizations of the noise ensemble. 



tracked again using Cauchy's argument principle (see Appendix [Cj . 



C. Arbitrary Couplings and Multiple Delays 

When there are multiple time delays involved in the synchronization or coordination process, in general, one cannot 
diagonalize the underlying system of coupled equations. This happens to be the case for the scenario with two 
types of time delay [Eq. ([28]) ] on unweighted (or globally weighted) graphs (as opposed to specific locally-weighted 
ones discussed in Sec. III.B). First, we briefly present a generally applicable method to determine the region of 
synchronizability/stability computationally [7y, [77]. For arbitrary couplings Cy , the deterministic part of Eq. (j2"5)l 
(from which one can extract the characteristic equation) becomes 

d t h t (t) = -dhi(t -r )+J2 C ijhj{t - r) , (39) 

3 

where Ci=^ ; Cu and r = r Q + r tr . After Laplace transform, these equations become 

shi(s) = -CMs)e- ST ° + Cijhj(s)e- ST , (40) 

3 

or equivalently, 

{sSij + CiSije-"-" ~ a je - ST ) hj(a) = . (41) 

3 

Hence, non-trivial solutions of the above system of equations require 

detM(s) = 0, (42) 

where 

My (a) = a5 tj + C ie ~ ST °5 l3 - C ZJ e~ ST . (43) 
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FIG. 12: (Color online) Comparison of (a) the widths and (b) the scaled widths for several coupling strengths a on a typical 
locally weighted BA network of size jV=100 and (k) ~ 6 for 7=0.2; simulated with D=l and At=0.001. 



(a) (b) 




FIG. 13: (Color online) Synchronization boundaries for several modes with (a) Xk < 1 and (b) > 1 of a weighted network, 
obeying Eq. (f3"5)) and determined by analyzing the zeros of Eq. (J37J). 



Stability or synchronizability requires that Re(s)<0 for all solutions of the above (transcendental) characteristic 
equation [Eq. (1421) ]. To identify the stability boundary of this coupled system, one does not need to know and determine 
the (infinitely many) complex solutions of the characteristic equation, but only whether all solutions have negative 
real parts. To test that, one again can employ the argument principle [76], (77| (Appendix [Cjl . Note that the above 
method can be immediately generalized to arbitrary heterogeneous (local and transmission) time delays. To compare 
synchronizability with locally weighted couplings of the same cost [Eq. here, we considered = uAij/ (k). The 
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FIG. 14: (Color online) Synchronization boundaries determined by analyzing the zeros of Eq. (|3T|) for various delay ratios 7, 
shown separately for (a) 7 < 0.5 and (b) 7 > 0.6. 




FIG. 15: (Color online) Width contributions for (a) several modes with 7 = 0.3 and (b) several delay ratios with A = 1.2, found 
by numerically integrating Eq. (135[) with D=l and a=l. The vertical lines correspond to the stability limits obtained from the 
analyses of the zeros of Eq. (|36p with the same A. 



results are shown in Fig. 1181 The synchronization boundary was determined using the above scheme, while the width 
was obtain by numerically integrating Eq. (1281) . Not only does local reweighting of the coupling strength improve 
synchronization, but it also extends the region of synchronizability. 
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FIG. 16: (Color online) Synchronization boundaries for typical (a) ER and (b) BA networks of several sizes with locally 
weighted couplings. The boundaries are found by numerical diagonalization and examining each mode through Eq. ()36[) . (c) 
Widths along a slice of constant r o =0.77 for the same iV=100 BA network used in (b). For stability comparison, the boundary 
is shown below with the slice indicated. 



IV. SUMMARY 



Through our investigations we have explored the impact and interplay of time delays, network structure, and 
coupling strength on synchronization and coordination in complex interconnected systems. Here, we considered only 
linear couplings, already yielding a rather rich phase diagram and response. While nonlinear effects are crucial in 
all real-life applications [9, I5ll - l53| , linearization and stability analysis about the synchronized state yields equations 
analogous to the ones considered here [H [1(3]. Hence, the detailed analysis of the linear problems can provide some 
insights to the complex phase diagrams and response of nonlinear problems. 

For a single uniform local delay, the synchronizability of a network is governed solely by the largest eigenvalue and 
the time delay. This result links the presence of larger hubs to the vulnerability of the system becoming unstable 
at smaller delays. The quality of synchronization within the stable regime is described by the width, which can be 
enumerated exactly for arbitrary symmetric couplings, provided the spectrum is known. We have also established 
the boundaries of the region of synhronizability in terms of the delay and the overall coupling strength (associated 
with communication rate) and provided the general scaling behavior of the width inside this region. Our results 
underscored the importance of the interplay of stochastic effects, network connections, and time delays, in that how 
"less" (in terms of local communication efforts) can be "more" efficient (in terms of global performance). 

For more general schemes with multiple time delays, we have shown how stability analysis in general delay differ- 
ential equations can be applied to ascertain the synchronizability of a network. For cases where, at least in principle, 
eigenmode decomposition is possible, we have identified the general scaling behavior of the width within the synchro- 
nizable regime. However, in these cases it is not always the same eigenvalue that determines stability for all 7. In 
the non-monotonic nature of the scaling function, we see that there is a fundamental limit to how well a network can 
synchronize in the presence of noise. In the case when transmission and reaction are two independent and significant 
sources of delay, there is an additional parameter for tuning: the ratio of local delay to the total delay. By fixing the 
local delay and cutting across different values of the ratio, there is the possibility that the network will enter into and 
emerge from synchronizable regions. Understanding these influences can guide network design in order to maintain 
and optimize synchronization by balancing the trade-offs in internodal communication and local processing. 
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FIG. 17: (Color online) The scaling functions of a typical locally weighted BA network and a typical ER network for two delay 
ratios, with both networks of size iV=100, found by numerically integrating Eq. (|34[) with D=l and At=0.001. The vertical 
lines correspond to the stability limits obtained from the analyses of the zeros of Eqs. (1371) . 
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Appendix A: Steady-state fluctuations for a single-variable stochastic delay equation 

For a single (linearized) stochastic variable h(t) with multiple time delays {r w }2 =1 and delta-correlated noise, one 
starts with the following general form 



d t h(t) = A h(t) + J2 A " h (t - t u ) + v(t) , 



(Al) 



where (rj(t)rj(t'))=2D&(t — t'). Formally, the noise r](t) plays the role of the inhomogeneous part of the above in- 
homogeneous linear first-order differential equation. Performing Laplace transformation [h(s) = / °° e~ st h(t)dt], the 
characteristic polynomial g(s) (and the corresponding characteristic equation) associated with the homogeneous (de- 
terministic) part of the above equation becomes 



g(s) = s - A - ^ A u e T " 



. 



(A2) 



For the initial condition h(t) =0 for t < (which we employ throughout this paper), we can easily obtain the Laplace 
transformed Green's function of Eq. (IA1[) (i.e., the solution when r)(t) is replaced by S(t — t')), which has the form 



G(s) 



e -st' 



(A3) 
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FIG. 18: (Color online) Scaled width curves for a typical BA network compared to those of a typical ER network of size 
N — 100 with (k) « 6 and D = 1, determined by numerically integrating Eq. (|28p for the two types of coupling schemes with 
7 = 0.1 and At = 0.01. 



Performing the inverse transform, one finds 



-i px +ioo -i pxo+ioo s (t-t') pS a {t—t') 



27r * Jxo-ioc 27Ti 7x Q - 4 oo 9(s) ^ g'(s a ) 

where s a (a — 1,2, . . .) are the zeros of the characteristic equation g(s) = on the complex plane [Eq. (| A2[) ] . In the 
above inverse transform, the infinite line of integration is parallel to the imaginary axis (s = xq) and is chosen to be 
to the right of all zeros of the characteristic polynomial in order to apply the residue theorem by closing the contour 
with an infinite semicircle to the left of this line. Note that the Green's function G(t, t') depends only on the variable 
t — t' , reflecting the time translation symmetry of the problem. Utilizing the Grenn's function, we can now formally 
write the general solution of Eq. (jAll) (for the same initial conditions) as 

h(t)= / dt'G(t,t')r,(t') = / dt'G(t,t')ri(t')= / dt' V ,, - r>(t') . (A5) 
Jo Jo Jo « 9 ( Sa > 

For more general initial conditions, see Ref. (78j . 

After averaging over the noise, one finds that the fluctuations of h(t) are 

«(*-*') ft ^ P s„(t-t") \ 



\Jo „ 9 0«) Jo p 9 M i 

ft ft ^_ e s a (t-t') e s e {t-t") 

Jo Jo „ 9 (s<x) V 9 (s/3) 



Y 1 f d f [ dt" e s ^ t ~ t '^e s ^ t - t "hD5(t' - t") 

^ 9 \s a )g {sf}) Jo Jo 



y 2D I dt > e (-+-/»)(t-f ) 

^ 9 (s a )g (S/9j Jo 



^ 9 \s a )g {sp){s a + sf}) 
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As is explicit from the above equation, the zero with the largest real part of all s a governs the long-time behavior of 
the stochastic variable h(t). In particular, h(t) reaches a stationary limit distribution for t — ¥ oo with a finite variance 
if and only if Re(s Q ) < for all a. In this case, 

(h 2 (oo)) = Y ,. . ,. . . , (A7) 

otherwise it diverges exponentially with time. Note that the condition for the existence of an asymptotic stationary 
limit distribution is the same as the one for the stability of the deterministic (homogeneous) part of Eq. (jAll) about 
the hi=0 fixed point UGH. 



Appendix B: Exact Scaling Functions for Time Delayed Stochastic Differential Equations 

Kiichler and Mensch [38| obtained the analytic stationary-state autocorrelation function for the stochastic delay- 
differential equation, 

d t h(t) = ah(t) + bh(t - r) + r)(t) , (Bl) 

with (r}(t)r}(t')) = 2D5{t — t'). Special cases (with suitably chosen coefficients a and b) can be directly utilized for two 
of our special cases in our network investigations. Specifically, for (i) unweighted networks with symmetric couplings 
and uniform local time delays with no transmission delays and in) the case with only transmission delays on the 
complete graph, to be discussed at the end of this Appendix, in Appendix IB 21 a and Appendix IB 21 b. respectively. 
Here, we briefly present an equivalent derivation of their results, using the formalism used in our paper. 
We define the stationary- state autocorrelation function as 

K(t) = (h{t')h(t' + t)) , (B2) 

where it is implicitly assumed that t' —¥ oo. From this definition and the invariance under time translation in the 
stationary state, it follows that the interpretation of the autocorrelation function can be formally extended to t<0 by 

K(t) = (h(t')h(t' + t)) = (h(t' + t)h(t')) = (h(t')h(t' - t)) = K(-t) , (B3) 

and also, 

k{t) = -k(-t) . (B4) 

As one would like to obtain a directly solvable equation of motion for the autocorrelation function, one must first 
find expressions for its time derivatives. Employing the equation of motion for h(t) [Eq. (|B1|) ]. we obtain for t > 
that 

K(t) = d t K{t) = d t (h(t')h(t' + t)) = (h(t')dth(t' + t)) = (h(t'){ah(t' + t) + bh{t' + t - t) + rj(t' + t)}) 

= a(h(t')h(t' + t)) + b{h{t')h{t' +t-r)) + (h(t')r]{t' + i)) (B5) 
= aK(t) + bK(t - t) , (B6) 

where in the last step we used (h(t')r)(t' + t)) =0 (i.e., Ito's convention [HD, [zl]). The above expression, combined 
with the (analytic) extension of the autocorrelation function in Eq. (|B3|) . yields the condition 

K(0) = aK(0) + bK{r) (B7) 



in the limit of t — > +0. Differentiating Eq. (|B6|) again with respect to t and exploiting the properties of Eqs. (|B3|) and 
(|B4|) . we find 

K(t) = aK(t) + bk(t -t) = aK(t) - bk(r - t) 

= a{aK(t) + bK(t - t)} - b{aK(r - t) + bK(-t)} 
= a{aK(t) + bK{t - t)} - b{aK(t - r) + bK(t)} 

= (a 2 - b 2 )K(t) . (B8) 
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Note that the reduction of the equation of motion of the autocorrelation function to a second order ordinary differential 
equation (with no delay) is a consequence of Eq. (|B1[) having only one delay time-scale. The general solution of Eq. (|B8[) 
can be written as 

K(t) = A cos(wi) + B sin(wi) (B9) 



with lo = \/b 2 — a 2 . From the definition of the autocorrelation function in Eq. (|B2[) and from some of the basic 
properties of the Green's function (see Appendix IB II below for details), it also follows [38| that 

K(0) = lim d t (h(t')h(t' + t)) = -D , (BIO) 

and from Eq. ([B7]). 

aK(0) + bK{r) = -D . (Bll) 

Thus, the second order ordinary differential equation Eq. (IB8|) with conditions Eqs. (|B10[) and (|B11[) can now be fully 
solved, yielding 

, , „ — lo + b sin(ojr) 

A = A' O = D— , B12 

Lu[a + bcos(ujT)} v ' 



K(0) D 

B = — — = . B13 

LO LO 



and 



Finally, the stationary-state variance of the stochastic variable governed by Eq. (|B1[) can be written as 

<*(«)) = (h(tMt)) = K(0) = D -" + bsi f T ] . (B14) 

LO[a + bcos{LOT)\ 

Following the aforementioned technical detours and details in Appendix IB II we will discuss in Appendix IB 21 the 
applications of the above result to obtain the scaling function of the fluctuations for the individual modes in specific 
networks. 



1. General properties of the autocorrelation function and the Green's Function 



From the definition of the autocorrelation function in Eq. (|B2j) and of the Green's function in Eq. (|A5j) . it follows 
that 

K[t) = (h(t')h(t' + *)) = ^jf duG(t',u)r)(u) + dv G{t' + t, «)»?(«) ^ 

{' rt'+t i-t' 
du dv G(t',u)G(t' + t,v)(r)(u)r)(v)) = 2D du G(t' ,u)G{t' + t,u) , (B15) 



and consequently 



K(t) = d t (h(t')h(t' + t)) =2Dd t [ duG{t',u)G{t' + t,u) = 2D [ du G(t' ,u)d t G{t' + t,u) 

Jo Jo 

= 2D [ duG(t' \u)(-d u )G(t' + t,u) = -2D [ du G(t' ,u)d u G(t' + t,u) . (B16) 
Jo Jo 

Hence, 

if(0) = -2D lim / duG(t',u)d u G{t' + t,u) = -2D [ du G(t' ,u)d u G(t' ,u) (B17) 
t ^°Jo Jo 

= -2D f du du * 3 ^'^ 2 = -D{G(t', t') - G(t', 0)} = -D{1 - 0} = -D , (B18) 

where in the second term of the last expression above we now explicitly exploited that G(t',t') = and G(t', 0) — > 
as t' — > oo. The former can be seen by a segment-by-segment integration and solution of Eq. (|B1[) with a delta source 
5(t — t') in the intervals (t — t') € [nr, (n + l)r], n = 0, 1,2, . . . [3g]; the solution in the [0, r] interval is particularly 
simple, G(t,t') — exp[a(t — £')]. The latter property is trivial in that the magnitude of the Green's function in the 
stationary state has to decay for large arguments. 
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2. Applications to Special Cases 



a. Unweighted Symmetric Couplings with Uniform Local Delays 

For symmetric couplings Cy with uniform local delays, the Laplacian I\y = &y J^i @il ~ Qj m Eq. (J7|) can, in 
principle, be diagonalized. Each mode is governed by Eq. ([5]), a special case of Eq. (|B1[) with a = 0, 6 = —A, and 
cj = |6| = A (A being the eigenvalue of the respective mode). From Eq. (|B14|) . the steady-state variance of each mode 
then reduces to 

<tf (oo)> = D±±^ = ^ri±^4 = DrfiXr) , (B19) 
Acos(At) Atcos(At) 

yielding the analytic scaling function for each mode 



1 + sin(aQ 

fix) = ft- , B20 

xcos(x) 



with the scaling variable x = At. 



b. Complete Graphs with Only Uniform Transmission Delays 



The exact stationary-state variance of Eq. (IB1|) can also be applied to complete graphs with global coupling er, 
which have no local delays but do have uniform transmission delays, i.e., Eq. (|30[) with 7 = 0, translating to a = —a, 
b = —cr/(N — 1) in Eq. (|B1I) . The analytic expression from Eq. (|B14[) for the stationary-state variance for each 
(non-uniform) mode becomes 

a + -jf—r sinh(ar) 

(h 2 (oo)) = D- ^-i -V-At , (B21) 

V V ;/ a[a+ 1 ^ T cosh(ar)] v ; 



with a = Va 2 - b 2 = a^Jl - 1/(N - l) 2 . 



Appendix C: Application of Cauchy's Argument Principle with Implementation 

For an arbitrary complex analytic function F(z), the number of zeros Nq inside a closed contour C (provided F(z) 
has no poles/singularities inside C) is given by Cauchy's argument principle (see, e.g., Ref. [80l|): 

N °=inlm dz= ^ AcargFizh (C1) 

where Ac engF(z) is the winding number of F(z) along the closed contour C. The characteristic equations studied 
in this paper can all be written as a sum of exponentials, hence there are no singularities. To determine the stability 
boundary, we follow Refs. [ifl and use Eq. (|C1|) to track the number of zeros of the characteristic equations with 
positive real part (i.e., on the positive real half plane) by substituting Eqs. (|5o| and P^l) for F(z). We employed a 
numerical algorithm [76| for enumerating the winding number with adaptive step size. We restate the method here: a 
step of size h along the contour in the direction i from s to s+hi is accepted if 9(s, s+(h/2)t < 9(s, s+hi) < e — 1 where 
8(s, s') = |arg(det M(s)) — arg(det M(s'))\ mod 2n. The subsequent step size is then h — > max{2, e/A}; unacceptable 
steps are retried with h — > h/2. The winding number is the count of the number of crossings of tt without a return in 
the opposite direction. 

We choose the contour so that it detects the first zero to cross the imaginary axis and acquire a positive real part. 
Note that the mode corresponding to the zero eigenvalue allows the solution z=0 for Eq. (|37|). so the zero at the origin 
is ignored. This can be easily achieved by choosing the left edge of the contour to be nonzero but still very small. This 
method can be applied to any network structure with any delay scheme, provided the approximate general behavior 
of the zeros is understood. 

As an example, consider the simplest system of two coupled nodes with uniform delay, which has a critical delay 
of 7r/4. Figure [in] shows two cases explored while finding the critical delay. In Fig. [TW ah t = ir/5 < t c and all real 
parts are non-positive so none fall within the contour. Tracking the argument (right column) shows that the winding 
number is correspondingly zero to verify that the delay is subcritical. Alternatively, r = 7r/3 > n c in Fig. I19f b) and 
there do indeed exist zeros with positive real parts that fall within the contour. The argument winds around the 
origin twice, signaling the presence of the first two zeros to cross the imaginary axis, indicating instability. 
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Re[s] 

FIG. 19: (Color online) Numerical integration ol Eq. (]C 1 1) to identity the presence ol zeros in the cases ol a system of two coupled 
nodes (r c =7r/4) for (a) r—n/5 and (b) r—n/3. The left column shows the zeros and the points sampled along the contour; the 
right column shows the argument of the characteristic function (angular coordinate) at these steps (radial coordinate). 

Appendix D: The Uniform Mode and the Width 

1. Eigenmode Decomposition 

In synchronization and coordination problems, it is natural to define an observable such as the width, which 
measures fluctuations with respect to the global mean, 

1 N 

^ 2 (t) = -J2Mt)-h(t)}\ (Dl) 

i=i 

where hit) = X^i^iW- -"- n wna t follows, we show that the amplitude associated with the uniform mode of the 
normalized Laplacian automatically drops out from the width. (In the case of unnormalized symmetric coupling, the 
expression for the width simplifies to the known form.) 

For our problem with two types of time delays and locally normalized couplings [Eq. (|34|) ]. decomposition along the 
right eigenvectors of K~ l L facilitates diagonalization. While this normalized Laplacian is a non-symmetric matrix, 
its eigenvalues are all real and non- negative (with the smallest being zero, Ao = 0). The corresponding (normalized) 
right eigenvector is 

|e )=iV- 1 /2(i,i i ...,if . (D2) 

Note that since the normalized Laplacian is non-symmetric, the eigenvectors are not orthogonal, i.e., (ei\ek) ^ $ik- 
To ease notational burden, in this subsection we use the bra-ket notation - not to be confused with ensemble average 
over the noise. In this notation, (-| is a row vector and |-) is a column vector, e.g., (eo| = 7V~ 1,/2 (1, 1, . . . , 1). Using 
this notation, the state vector is denoted by 

\h(t)) = (h 1 (t),h 2 (t),...,h N (t)) T , (D3) 

while the state vector relative to the mean is 

\h(t)-h(t)) = (hi(t) -h(t),h 2 (t) -h(t),...,h N (t) -h{t)) T 
= {hi(t), h 2 (t), h N (t)) T - h(t)(l, 1, . . . , 1) T 

= \h(t))-h(t)VN\e ) = (l-\e )(e \)\h{t)) . (D4) 
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Employing the above formalism, the width can be written as 

N 



i=l 

Now we express the state vector as the linear combination of the eigenvectors of the underlying Laplacian, 

JV-l 



\h(t)) J2 hk(t)\e k ) . (D6) 



fc=0 

Employing the above eigenmode decomposition, (h — h\h — h) can be written as 

N-l N-l 



(h - h\h - h) = (h\(l - \e )(e \) 2 \h) = (h\(l - \e )(e \)\h) = £ h k (t){e k \ (1 - |e )(e |) £ fc,(t)|e,) 

k=a 1=0 

N-l N-l 

= ^2 hk{t)hi{t)(e k \ (1 - |eo)<eo|) |ej) = ^ h{t)hi(t) ((e fe |ej) - (e fc |e )(e |e ; )) 



fc,Z=0 fc,Z=0 

- 5] h k {t)hi{t) ((e k \e t ) - (e fc |e )(e |e ; )) = ^ h k {t)hi{t) (E kl - E ko E ol ) , (D7) 

where Ejy = (e/-|e;). As can be seen explicitly from Eq. (|D7|) . the terms where either A: or I are zero drop out from 
the sum (as Eqq=1). It is also clear from Eq. (|D7|) that (/i — h\h — h) = J2 k ^o wnen the underlying coupling is 
symmetric (and consequently the eigenvectors form an orthogonal set, E k i=8 k i). Finally, the width can be written as 

1 N 1 1 

w 2 (t) = - X>iW - W = ^ - % - h) = - J2 hk(t)hi(t) (E kl - E k0 E i) . (D8) 

Note that the above result can be immediately applied to the case of symmetric coupling with no transmission delays 
[Eq. ||7J]. There, the eigenvectors of the corresponding Laplacian form an orthogonal set, and the above expression 

collapses to w 2 {t) = ± J2 k =i H (*) 0- 

2. Ensemble Average over the Noise 

We now use the general form of the solution given in Appendix [A] [Eq. (|A5|) ] for the respective eigenmodes of 
normalized Laplacian coupling with two types of time delays [Eq. (|35[) ] , giving 

r l _ e s ka (t-t') 

h k (t) = / dt'Y, -, > . (D9) 

where Sfe Q is the ath solution of the kth mode for the characteristic equation g k (s) = [Eq. After averaging 

over the noise, one obtains for the two-point function 

„ (I — e (Sfea+S W )f) 

(h k (t)h(t)) = -2D X H £ v S v w ■ (D10) 

^ 9 k {Ska)g l (Sll3){S ka + Sip) 

In the stationary state, one must have Re(sfc Q ) < for all k and a. Thus, the stationary state width can be written 
as 

<- 2 M> = lim 1 2 <M*)M*)> - EM .^EE ( g »7^)x« . (D11) 
t ^°° N k % N k %t? 9 k (sk a ) gi ( S i(3)(s ka + sip) 
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